---
title: 'Replication File for "Social Disruption, Gun-Buying, and Anti-System Beliefs" '
author: "Matthew J. Lacombe, Matthew D. Simonson, Jon Green, James N. Druckman"
date: "`r Sys.Date()`"
output: 
  html_document:
    toc: true
    toc_depth: 2
---

# Introduction
Welcome to the replication file! This file draws on data from `lacombe_disruption_data.rds` in the same directory as this file and outputs to subdirectories named `figures`, `tables`, and `appendix_tables`.  This script was created in RStudio 2022.07.1 running R 4.1.1. Please contact [Matthew Simonson](mailto:msimonson@sas.upenn.edu) with questions.

# Setup

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
cat("\14")

dir.create(path = "tables")
dir.create(path = "appendix_tables")
dir.create(path = "figures")
```

## Libraries

```{r, message=F}
# Standard
library(tidyverse)
library(data.table)

# For Regressions
library(fixest)
library(broom)

# Four Surveys
library(survey)
library(srvyr)

# For Plots
library(viridis)
library(ggforestplot)

# For Tables
library(kableExtra)
library(xtable)
library(gt)
library(gtsummary)
library(modelsummary)
```

## Data
```{r}
d <- readRDS("lacombe_disruption_data.rds")


s <- d %>% as_survey_design(weight = weights)

# Subsets
d.no.new.owners <- d %>% filter(new_owner == F)
d.all.owners <- d %>% filter(gun_owner == T)
d.new.and.oldnobuy <- d %>% filter(gun_owner == T, new_owner == T | gun_pan == F)
d.buyers <- d %>% filter(gun_pan == T)
```

## Functions

```{r}
tidy_one_predictor <- function(model, dv, predictor, ci.level = 0.95, 
                               inner.ci.level = 0.9, double.ci = T){
  
  output <- tidy(model, conf.int = T, conf.level = ci.level) %>% 
    filter(term == predictor) %>% 
    mutate(dv = dv) %>% relocate(dv) %>% 
    rename("predictor" = term)
  
  if(double.ci == T){
    inner.ci <- tidy(model, conf.int = T, conf.level = inner.ci.level) %>% 
      filter(term == predictor) %>% 
      transmute(conf.low.inner = conf.low, conf.high.inner = conf.high)
    output <- cbind(output, inner.ci)
  }
  output
}

tidy_many_outcomes <- function(model_list, outcome_list, predictor, 
                               ci.level = 0.95, inner.ci.level = 0.9, 
                               double.ci = T){
  map2(model_list, outcome_list, tidy_one_predictor, predictor = predictor, 
       ci.level = ci.level, inner.ci.level = inner.ci.level, double.ci = double.ci) %>% 
    bind_rows()
}

tidy_many_predictors <- function(model, dv, predictor_list, ci.level = 0.95, 
                                 inner.ci.level = 0.9, double.ci = T){
  
  output <- tidy(model, conf.int = T, conf.level = ci.level) %>% 
    filter(term %in% predictor_list) %>% 
    mutate(dv = dv) %>% relocate(dv) %>% 
    rename("predictor" = term)
  
  if(double.ci == T){
    inner.ci <- tidy(model, conf.int = T, conf.level = inner.ci.level) %>% 
      filter(term %in% predictor_list) %>% 
      transmute(conf.low.inner = conf.low, conf.high.inner = conf.high)
    output <- cbind(output, inner.ci)
  }
  output
}
```

## Fixest Settings
```{r}

setFixest_etable(digits = 2)

# set the control variables for feglm regressions
setFixest_fml(..controls = ~
                Black + Asian_Am + Hispanic + Other_Race +
                Female + Parent +
                age_in_decades +
                college +
                income_in_10ks +
                Rural + Suburban +
                white_evang +
                Democrat + Republican + Other_Party +
                ideology +
                region_8 + 
                ongoing_hardship_index + covid_house)

# set names for regression tables
setFixest_dict(c(gun_pan = "Pandemic Gun-Buyer",
                existing_owner = "Pre-Existing Gun Owner",
                gun_owner = "Gun Owner",
                new_owner = "New Gun Owner",
                
                trump_win = "Trump Won",
                vaccine_misinfo_index = "Vaccine Conspiracy Index",
                vac_minfo_alter_dna = "Alters DNA",
                vac_minfo_fetal_tissue = "Contains Fetal Tissue",
                vac_minfo_microchips = "Contains Microchips",
                vac_minfo_infertility = "Causes Infertility",
                vac_minfo_clinically_tested_NO = "Not Tested",
                trust_news_media = "Trust News Media",
                trust_gov_health_index = "Trust Health Officials",
                trust_scientists = "Trust Scientists",
                reason_hobby = "Hobbyist Reasons",
                reason_threat = "Threat Reasons",
                reason_other = "Other",
              
                "Black" = "Black", 
                 "Asian_Am" = "Asian",
                 "Hispanic" = "Hispanic" , 
                 "Other_Race" = "Other Race" ,
                 "Female" = "Female" , 
                 "Parent" = "Children in HH" ,
                 "age_in_decades" = "Age (Decades)" ,
                 "college" = "College" ,
                 "income_in_10ks" = "HH Income (10k)" ,
                 "Rural" = "Rural" ,
                 "Suburban" = "Suburban" ,
                 "white_evang" = "White Evangelical" ,
                 "Democrat" = "Democrat" , 
                 "Republican" = "Republican" , 
                 "Other_Party" = "Other Party" ,
                 "ideology" = "Ideological Identity" ,
                 "region_8Rockies" = "Region: Rockies" , 
                 "region_8Southwest" = "Region: Southwest",
                 "region_8Great Plains" = "Region: Great Plains",
                 "region_8South" = "Region: South",
                 "region_8Midwest" = "Region: Midwest",
                 "region_8Mid-Atlantic" = "Region: Mid-Atlantic",
                 "region_8New England" = "Region: New England",
                 "ongoing_hardship_index" = "Hardship Index" , 
                 "covid_house" = "COVID in HH"))
cmap <- c("Black" = "Black", 
          "Asian_Am" = "Asian",
          "Hispanic" = "Hispanic" , 
          "Other_Race" = "Other Race" ,
          "Female" = "Female" , 
          "Parent" = "Children in HH" ,
          "age_in_decades" = "Age (Decades)" ,
          "college" = "College" ,
          "income_in_10ks" = "HH Income (10k)" ,
          "Rural" = "Rural" ,
          "Suburban" = "Suburban" ,
          "white_evang" = "White Evangelical" ,
          "Democrat" = "Democrat" , 
          "Republican" = "Republican" , 
          "Other_Party" = "Other Party" ,
          "ideology" = "Ideological Identity" ,
          "region_8Rockies" = "Region: Rockies" , 
          "region_8Southwest" = "Region: Southwest",
          "region_8Great Plains" = "Region: Great Plains",
          "region_8South" = "Region: South",
          "region_8Midwest" = "Region: Midwest",
          "region_8Mid-Atlantic" = "Region: Mid-Atlantic",
          "region_8New England" = "Region: New England",
          "ongoing_hardship_index" = "Hardship Index" , 
          "covid_house" = "COVID in HH")
```

# Who Bought Guns

Regressions
```{r}
mod.who.bought <- feglm(gun_pan ~ existing_owner + ..controls, 
                       data = d, family = gaussian(), se = "hetero")
mod.who.bought.bin <- feglm(gun_pan ~ existing_owner + ..controls, 
                        data = d, family = "binomial", se = "hetero")
```

Table C1
```{r}
etable("OLS" = mod.who.bought, 
       "Logit" = mod.who.bought.bin, 
       adjustbox = .45,
       replace = T, 
       title = "Bought Guns During the Pandemic",
       file = "appendix_tables/c1.tex")
```

## Figure 1
```{r}
df.who.bought <- tidy_many_predictors(model = mod.who.bought, dv = "gun_pan",
                     predictor_list = list(
                       "Black", "Asian_Am", "Hispanic", #"Other_Race",
                       "Female", "Parent",
                       "age_in_decades",
                       "college",
                       "income_in_10ks",
                       "Rural", "Suburban",
                       "white_evang",
                       "Democrat", "Republican",# "Other_Party",
                       "ideology",
                       "region_8", 
                       "ongoing_hardship_index", "covid_house"
                     )) %>%
  mutate(significant = 1*(p.value < 0.05),
         predictor = fct_recode(
           predictor, 
           "Asian" = "Asian_Am",
           #"Other Race/Ethnicty" = "Other_Race",
           "Age (in decades)" = "age_in_decades",
           "College" = "college",
           "Income (in $10,000)" = "income_in_10ks",
           "Rural (vs. Urban)" = "Rural",
           "Suburban (vs. Urban)" = "Suburban",
           "White Evangelical" = "white_evang",
           "Democrat (vs. Independent)" = "Democrat",
           "Republican (vs. Independent)" = "Republican",
           #"Other Party (vs. Independent)" = "Other_Party",
           "Ideology" = "ideology",
           "Hardship Index" = "ongoing_hardship_index",
           "House with Covid-19" = "covid_house") %>% fct_reorder(estimate),
         threat = if_else(predictor %in% c("Hardship Index", "House with Covid-19"), 1, 0))

figure_1 <- 
  ggplot(data = df.who.bought, 
         aes(y=predictor, color = threat)) + 
  geom_pointrange(aes(x = estimate, 
                      xmin = conf.low, 
                      xmax = conf.high)) +
  geom_linerange(aes(xmin = conf.low.inner, 
                     xmax = conf.high.inner), 
                 size = 1.5) +
  geom_vline(xintercept = 0, lty = "dashed") + 
  theme_forest() + 
  geom_stripes(odd = "#22222222", 
               even = "#00000000") + 
  labs(x = "Estimate", y = "Variable", 
       title = "Figure 1: Marginal Probability of Purchasing a Gun during the Pandemic",
       subtitle = "Estimates for threat variables (hardship index and COVID-19 in household) highlighted",
       caption = "Not shown: intercept, region, pre-existing ownership status, other party, other race/ethnicity") +
  scale_color_gradient(low = "grey", 
                       high = "black") +
  theme(legend.text = 
          element_text(size=12), 
        legend.position = "none")

# All the plot sizes in this script may need adjusting
ggsave(figure_1, file = "figures/Figure_1.pdf", width = 10, height = 5)

```

# Conspiracy Beliefs

## Subset excluding new owners
```{r}
left.conspiracy <- lapply(c("trump_win", "vaccine_misinfo_index",
                         "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                         "vac_minfo_microchips", "vac_minfo_infertility", 
                         "vac_minfo_clinically_tested_NO"), function(x){
                           fixest::feglm(as.formula(paste0(x, "~ existing_owner + ..controls")),
                                         data = d.no.new.owners, 
                                         family = gaussian(), 
                                         se = "hetero")
                         })

left.conspiracy.bin <- lapply(c("vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                         "vac_minfo_microchips", "vac_minfo_infertility", 
                         "vac_minfo_clinically_tested_NO"), function(x){
                           fixest::feglm(as.formula(paste0(x, "~ existing_owner + ..controls")),
                                         data = d.no.new.owners, 
                                         family = "binomial", 
                                         se = "hetero")
                         })
```

Tables C2 and C3
```{r}
etable(left.conspiracy,
       title = "Conspiracy Beliefs by Gun Ownership (Pre-Existing Owners vs. Non-Owners) (OLS)",
       adjustbox = T,
       replace = T,
       file = "appendix_tables/c2.tex")
etable(left.conspiracy.bin,
       title = "Conspiracy Beliefs by Gun Ownership (Pre-Existing Owners vs. Non-Owners) (Logit)",
       adjustbox = T,
       replace = T,
       file = "appendix_tables/c3.tex")
```

## Subset of all owners
```{r}
left.conspiracy.df <- tidy_many_outcomes(left.conspiracy, 
                                      list("trump_win", "vaccine_misinfo_index",
                                           "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                                           "vac_minfo_microchips", "vac_minfo_infertility", 
                                           "vac_minfo_clinically_tested_NO"),
                                      "existing_owner") %>%
  mutate(heading = "Effect of Pre-Existing Gun Owner \n(vs. Non-Owner)")


right.conspiracy <- lapply(c("trump_win", "vaccine_misinfo_index",
                          "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                          "vac_minfo_microchips", "vac_minfo_infertility", 
                          "vac_minfo_clinically_tested_NO"), function(x){
                            fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                          data = d.all.owners, 
                                          family = gaussian(), 
                                          se = "hetero")
                          })


right.conspiracy.bin <- lapply(c("vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                             "vac_minfo_microchips", "vac_minfo_infertility", 
                             "vac_minfo_clinically_tested_NO"), function(x){
                               fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                             data = d.all.owners, 
                                             family = "binomial", 
                                             se = "hetero")
                             })
```

Tables C4 and C5
```{r}

etable(right.conspiracy,
       title = "Conspiracy Beliefs by Gun Ownership (New Owners vs. Pre-Existing Owners) (OLS)",
       adjustbox = T,
       replace = T,
       file = "appendix_tables/c4.tex")

etable(right.conspiracy.bin,
       title = "Conspiracy Beliefs by Gun Ownership (New Owners vs. Pre-Existing Owners) (Logit)",
       adjustbox = T,
       replace = T,
       file = "appendix_tables/c5.tex")
```

## Subset of owners excluding pandemic buyers
```{r}
right.conspiracy.df <- tidy_many_outcomes(right.conspiracy, 
                                       list("trump_win", "vaccine_misinfo_index",
                                            "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                                            "vac_minfo_microchips", "vac_minfo_infertility", 
                                            "vac_minfo_clinically_tested_NO"),
                                       "new_owner") %>% 
  mutate(heading = "Effect of New Gun Owner \n(vs. Pre-Existing Gun Owner)")

# check pval for trump_win
summary(right.conspiracy[[1]])

farright.conspiracy <- lapply(c("trump_win", "vaccine_misinfo_index",
                          "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                          "vac_minfo_microchips", "vac_minfo_infertility", 
                          "vac_minfo_clinically_tested_NO"), function(x){
                            fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                          data = d.new.and.oldnobuy, 
                                          family = gaussian(), 
                                          se = "hetero")
                          })
farright.conspiracy.bin <- lapply(c("vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                              "vac_minfo_microchips", "vac_minfo_infertility", 
                              "vac_minfo_clinically_tested_NO"), function(x){
                                fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                              data = d.new.and.oldnobuy, 
                                              family = "binomial", 
                                              se = "hetero")
                              })
# regression tables
```

Tables C6 and C7
```{r}

etable(farright.conspiracy,                           
       title = "Conspiracy Beliefs by Gun Ownership (New Gun Owners vs. Pre-Existing Gun Owners With No Pandemic Purchases) (OLS)",
       adjustbox = T,
       replace = T,
       file = "appendix_tables/c6.tex")

etable(farright.conspiracy.bin,  
       title = "Conspiracy Beliefs by Gun Ownership (New Gun Owners vs. Pre-Existing Gun Owners With No Pandemic Purchases) (Logit)",
       adjustbox = T,
       replace = T,
       file ="appendix_tables/c7.tex")
```

```{r}
farright.conspiracy.df <- tidy_many_outcomes(
  farright.conspiracy, 
  list("trump_win", "vaccine_misinfo_index",
       "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
       "vac_minfo_microchips", "vac_minfo_infertility", 
       "vac_minfo_clinically_tested_NO"),
  "new_owner") %>% 
  mutate(heading = "Effect of New Gun Owner \n(vs. Pre-Existing Gun Owner who \ndid not buy in pandemic)")
```

## Figure 2
```{r}
conspiracy.df <- bind_rows(left.conspiracy.df, 
                        right.conspiracy.df, 
                        farright.conspiracy.df) %>%
  mutate(dv = fct_rev(fct_inorder(fct_recode(
    dv, 
    "Trump Won 2020 Election" = "trump_win", 
    "Vaccine Conspiracy Index" = "vaccine_misinfo_index",
    "Conspiracy: Vaccine Alters DNA" = "vac_minfo_alter_dna",
    "Conspiracy: Vaccine Contains Fetal Tissue" = "vac_minfo_fetal_tissue",
    "Conspiracy: Vaccine Contains Microchips" = "vac_minfo_microchips",
    "Conspiracy: Vaccine Causes Infertility" = "vac_minfo_infertility",
    "Conspiracy: Vaccine Not Tested" = "vac_minfo_clinically_tested_NO"))),
    heading = fct_inorder(heading))

  
figure_2 <- ggplot(data = conspiracy.df) + 
  geom_pointrange(aes(y=dv, x = estimate, 
                      xmin = conf.low, xmax = conf.high)) +
  geom_linerange(aes(y=dv, 
                     xmin = conf.low.inner, 
                     xmax = conf.high.inner), 
                 size = 1.5) +
  geom_vline(xintercept = 0, lty = "dashed") + 
  facet_wrap(vars(heading)) +  
  theme_forest() + 
  labs(x = "Estimate", y = "Outcome",
       title = "Figure 2: Conspiracy beliefs among existing and new gun owners",
       subtitle = "Estimates from linear regressions with controls (not shown). 90% and 95% uncertainty intervals use robust standard errors.")

ggsave(figure_2, file = "figures/Figure_2.pdf", width = 13, height = 5)

```

# Trust

## Subset excluding new owners
```{r}
left.trust <- lapply(c("trust_news_media", "trust_gov_health_index",
                        "trust_scientists"), function(x){
                               fixest::feglm(as.formula(paste0(x, "~ existing_owner + ..controls")),
                                             data = d.no.new.owners, 
                                             family = gaussian(), 
                                             se = "hetero")
                             })
```

Table C8
```{r}

etable(left.trust, 
       adjustbox = .7,
       replace = T,
       title = "Trust by Gun Ownership (Pre-Existing Owners vs. Non-Owners)",
       file = "appendix_tables/c8.tex")
```


```{r}
left.trust.df <- tidy_many_outcomes(left.trust, 
                                      list("trust_news_media", "trust_gov_health_index", "trust_scientists"),
                                      "existing_owner") %>%
  mutate(heading = "Effect of Pre-Existing Gun Owner \n(vs. Non-Owner)")
```

##Subset of all owners
```{r}
right.trust <- lapply(c("trust_news_media", "trust_gov_health_index",
                       "trust_scientists"), function(x){
                         fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                       data = d.all.owners, 
                                       family = gaussian(), 
                                       se = "hetero")
                       })
```

Table C9
```{r}

etable(right.trust,
       adjustbox = .7,
       replace = T,
       title = "Trust by Gun Ownership (New Owners vs. Pre-Existing Owners)",
       file = "appendix_tables/c9.tex")
```


```{r}
right.trust.df <- tidy_many_outcomes(right.trust, 
                                       list("trust_news_media", "trust_gov_health_index", "trust_scientists"),
                                       "new_owner") %>% 
  mutate(heading = "Effect of New Gun Owner \n(vs. Pre-Existing Gun Owner)")
```

## Subset of owners excluding pandemic buyers
```{r}
farright.trust <- lapply(c("trust_news_media", "trust_gov_health_index",
                        "trust_scientists"), function(x){
                          fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                        data = d.new.and.oldnobuy, 
                                        family = gaussian(), 
                                        se = "hetero")
                        })
```

Table C10
```{r}
etable(farright.trust, 
       title = "Trust by Gun Ownership (New Owners vs. Pre-Existing Owners With No Pandemic Purchases)", 
       adjustbox = .7,
       replace = T,
       file = "appendix_tables/c10.tex")
```


```{r}
farright.trust.df <- tidy_many_outcomes(farright.trust, 
                                     list("trust_news_media", "trust_gov_health_index", "trust_scientists"),
                                     "new_owner") %>% 
  mutate(heading = "Effect of New Gun Owner \n(vs. Pre-Existing Gun Owner who \ndid not buy in pandemic)")
```

## Figure 3

```{r}
trust.df <- bind_rows(left.trust.df, right.trust.df, farright.trust.df) %>%
  mutate(dv = fct_rev(fct_inorder(fct_recode(
    dv, 
    "Trust Health Officials" = "trust_gov_health_index", 
    "Trust Scientists" = "trust_scientists",
    "Trust News Media" = "trust_news_media"))),
    heading = fct_inorder(heading))

figure_3 <- 
  ggplot(data = trust.df) + 
  geom_pointrange(aes(y=dv, 
                      x = estimate, 
                      xmin = conf.low, 
                      xmax = conf.high)) +
  geom_linerange(aes(y=dv, 
                     xmin = conf.low.inner, 
                     xmax = conf.high.inner), 
                 size = 1.5) +
  geom_vline(xintercept = 0, lty = "dashed") + 
  facet_wrap(vars(heading)) +  
  theme_forest()  + 
  labs(x = "Estimate", y = "Outcome",
       title = "Figure 3: Institutional trust among existing and new gun owners",
       subtitle = "Linear regression with controls (not shown). 90% and 95% uncertainty intervals use robust standard errors.")

ggsave(figure_3, file = "figures/Figure_3.pdf", width = 14, height = 5)

```

# Reasons

## Aggregated Reasons
```{r}
mod.reasons.agg <- lapply(c("(reason_hunt == 1 | reason_shoot == 1)", 
                        "(reason_crime == 1|reason_gov == 1|reason_covid == 1 |reason_lockdown == 1|reason_election == 1|reason_protection == 1)",
                        "reason_other"), function(x){
                          fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                        data = d.buyers, 
                                        family = gaussian(), 
                                        se = "hetero")
                        })

mod.reasons.agg.bin <- lapply(c("(reason_hunt == 1 | reason_shoot == 1)", 
                            "(reason_crime == 1|reason_gov == 1|reason_covid == 1 |reason_lockdown == 1|reason_election == 1|reason_protection == 1)",
                            "reason_other"), function(x){
                              fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                            data = d.buyers, 
                                            family = "binomial", 
                                            se = "hetero")
                            })
```

Tables C11 and C12
```{r}
etable(mod.reasons.agg, 
       adjustbox = .4,   
       replace = T,
       depvar = F,
       headers = c("Hobby", "Threat", "Other"),
       title = "Aggregated Reasons for Pandemic Gun Purchases (New Owners vs. Pre-Existing Owners) (OLS)",
       file = "appendix_tables/c11.tex")

etable(mod.reasons.agg.bin, 
       adjustbox = .4,
       replace = T,
       depvar = F,
       headers = c("Hobby", "Threat", "Other"),
       title = "Aggregated Reasons for Pandemic Gun Purchases (New Owners vs. Pre-Existing Owners) (Logit)",
       file = "appendix_tables/c12.tex")
```

## Figure 4
```{r}
mod.reasons.agg.df <- tidy_many_outcomes(mod.reasons.agg, 
                                     list("(reason_hunt == 1 | reason_shoot == 1)",
                                          "(reason_crime == 1 | reason_gov == 1 | reason_covid == 1 | reason_lockdown == 1 | reason_election == 1 | reason_protection == 1)",
                                          "reason_other"),
                                     predictor = "new_owner") %>%
  mutate(dv = fct_rev(fct_inorder(fct_recode(dv, "Hobbyist Reason(s)" = "(reason_hunt == 1 | reason_shoot == 1)", 
                                             "Threat Reason(s)" = "(reason_crime == 1 | reason_gov == 1 | reason_covid == 1 | reason_lockdown == 1 | reason_election == 1 | reason_protection == 1)", 
                                             "Other" = "reason_other"))),
         header = "Effect of New Gun Owner (vs. Pre-Existing Gun Owner)
Sample: Respondents who bought guns from March 2020 onwards")

figure_4 <- 
  ggplot(data = mod.reasons.agg.df) + 
  geom_pointrange(aes(y=dv, 
                      x = estimate, 
                      xmin = conf.low, 
                      xmax = conf.high)) +
  geom_linerange(aes(y=dv, 
                     xmin = conf.low.inner, 
                     xmax = conf.high.inner), 
                 size = 1.5) +
  geom_vline(xintercept = 0, lty = "dashed") + 
  theme_forest()  +
  facet_wrap(vars(header))+
  scale_x_continuous(name = "Estimate",
                     breaks = c(-0.3, -0.25, -0.2, -0.15, -0.1, -0.05,
                                0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3),
                     limits = c(-0.35, 0.35),
                     labels = c("-0.3","-0.25", "-0.2",
                                "0.15\n\nLikelier to be selected by pre-existing owners",
                                "-0.1", "0.05", "0", "0.05", "0.1",
                                "0.15\n\nLikelier to be selected by new owners",
                                "0.2","0.25", "0.3"))+
  labs(title = "Figure 4: Reasons for purchasing guns during the pandemic",
       subtitle = "Estimates from linear regression with controls (not shown). 90% and 95% uncertainty intervals use robust standard errors.",
       y = "Outcome",
       caption = 'Hobbyist reasons: "hunting", "target shooting"\nThreat reasons: "protection against crime," "protection against the government," "because of COVID-19," "because of lockdown and restrictions," "because of the election," "protection against someone I know personally"')+
  theme(plot.caption = element_text(hjust = 0, size = 6))

ggsave(figure_4, file = "figures/Figure_4.pdf", width = 12, height = 5)

```

## Alphas
```{r}

alpha.conspiracy <- 
  as_tibble(d) %>%
  dplyr::select(vac_minfo_alter_dna, 
                vac_minfo_fetal_tissue, 
                vac_minfo_microchips, 
                vac_minfo_infertility, 
                vac_minfo_clinically_tested_NO) %>%
  mutate_all(.funs = as.numeric) %>%
  psych::alpha()

alpha.health.inst <- 
  as_tibble(d) %>%
  dplyr::select(trust_CDC, 
                trust_fauci, 
                trust_FDA) %>%
  mutate_all(.funs = as.numeric) %>%
  psych::alpha()

alpha.hardship <- 
  as_tibble(d) %>%
  dplyr::select(starts_with("ongoing")) %>%
  dplyr::select(-ends_with("_index")) %>%
  psych::alpha()

alpha.threat.reason <-
  as_tibble(d) %>%
  dplyr::select(reason_crime, 
                reason_gov, 
                reason_covid, 
                reason_lockdown, 
                reason_election, 
                reason_protection) %>%
  filter(!is.na(reason_crime)) %>%
  psych::alpha()

alpha.hobby.reason <-
  as_tibble(d) %>%
  dplyr::select(reason_hunt,
                reason_shoot) %>%
  filter(!is.na(reason_hunt)) %>%
  psych::alpha()
```

# Robustness Checks

## Table 1: Conspiracy Beliefs
```{r}
stars <- c("***"=0.01, "**"=0.05, "*"=0.10)

robustness.misinfo <- lapply(c("trump_win", "vaccine_misinfo_index",
                         "vac_minfo_alter_dna", "vac_minfo_fetal_tissue", 
                         "vac_minfo_microchips", "vac_minfo_infertility", 
                         "vac_minfo_clinically_tested_NO"), function(x){
                           fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                         data = d.all.owners, 
                                         family = gaussian(), 
                                         se = "hetero")
                         })


names.rob.misinfo <- c("Trump Won","Vaccine Conspiracy Index")
names(robustness.misinfo)[c(1,2)] <- names.rob.misinfo

modelsummary::modelsummary(robustness.misinfo[c(1,2)],
                           coef_map = cmap,
                           gof_map = c("nobs","pseudo.r.squared",
                                       "AIC","BIC","logLik",
                                       "vcov.type"),
                           stars = c("*" = .1, "**" = .05, "***" = .001),
                           title = "Correlates of Conspiracy Beliefs (All Owners)",
                           output = "tables/table_1.html")

etable(robustness.misinfo[c(1,2)], 
       title = "Correlates of Conspiracy Beliefs (All Owners)",
       replace = T,
       adjustbox = .5,
       file = "appendix_tables/c13.tex")
```

## Table 2: Trust
```{r}
robustness.trust <- lapply(c("trust_news_media", "trust_gov_health_index",
                       "trust_scientists"), function(x){
                         fixest::feglm(as.formula(paste0(x, "~ new_owner + ..controls")),
                                       data = d.all.owners, 
                                       family = gaussian(), 
                                       se = "hetero")
                       })

names.rob.trust<- paste("Trust:",c("News Media","Health Index","Scientists"))
names(robustness.trust) <- names.rob.trust

modelsummary::modelsummary(robustness.trust,
                           coef_map = cmap,
                           gof_map = c("nobs","pseudo.r.squared",
                                       "AIC","BIC","logLik",
                                       "vcov.type"),
                           stars = c("*" = .1, "**" = .05, "***" = .001),
                           title = "Correlates of Institutional Trust (All Owners)",
                           output = "tables/table_2.html")

etable(robustness.trust, 
       title = "Correlates of Institutional Trust (All Owners)",
       replace = T,
       adjustbox = .5,
       file = "appendix_tables/c14.tex")
```

# Summary Statistics

## Descriptives Table

```{r}
out1 <- tbl_svysummary(s, 
                       statistic = list(all_continuous() ~ "{mean}", 
                                        all_categorical() ~ "{p}%"),
                       include = c("gender", "race", "urban_type", "age", "edu_", "income_brackets_4", "party", "ideo_fct_3", "trump_won", "existing_owner",
                                   "new_owner"),
                       label = list(
                         edu_ ~ "Education",
                         urban_type ~ "Community",
                         income_brackets_4 ~ "Income Bracket",
                         ideo_fct_3 ~ "Ideology",
                         trump_won ~ "Trump won",
                         party ~ "Party",
                         age ~ "Age",
                         race ~ "Race",
                         gender ~ "Gender",
                         existing_owner ~ "Pre-Existing Owner",
                         new_owner ~ "New Owner"
                       ),
                       missing = "no")

out2 <- tbl_svysummary(s %>% filter(existing_owner == 1), 
                       statistic = list(all_continuous() ~ "{mean}", 
                                        all_categorical() ~ "{p}%"),
                       include = c("gender", "race", "urban_type", "age", "edu_", "income_brackets_4", "party", "ideo_fct_3", "trump_won", "existing_owner",
                                   "new_owner"),
                       label = list(
                         edu_ ~ "Education",
                         urban_type ~ "Community",
                         income_brackets_4 ~ "Income Bracket",
                         ideo_fct_3 ~ "Ideology",
                         trump_won ~ "Trump won",
                         party ~ "Party",
                         age ~ "Age",
                         race ~ "Race",
                         gender ~ "Gender",
                         existing_owner ~ "Pre-Existing Owner",
                         new_owner ~ "New Owner"
                       ),
                       missing = "no")

out3 <- tbl_svysummary(s %>% filter(gun_pan == 1), 
                       statistic = list(all_continuous() ~ "{mean}", 
                                        all_categorical() ~ "{p}%"),
                       include = c("gender", "race", "urban_type", "age", "edu_", "income_brackets_4", "party", "ideo_fct_3", "trump_won", "existing_owner",
                                   "new_owner"),
                       label = list(
                         edu_ ~ "Education",
                         urban_type ~ "Community",
                         income_brackets_4 ~ "Income Bracket",
                         ideo_fct_3 ~ "Ideology",
                         trump_won ~ "Trump won",
                         party ~ "Party",
                         age ~ "Age",
                         race ~ "Race",
                         gender ~ "Gender",
                         existing_owner ~ "Pre-Existing Owner",
                         new_owner ~ "New Owner"
                       ),
                       missing = "no")

out4 <- tbl_svysummary(s %>% filter(new_owner == 1), 
                       statistic = list(all_continuous() ~ "{mean}", 
                                        all_categorical() ~ "{p}%"),
                       include = c("gender", "race", "urban_type", "age", "edu_", "income_brackets_4", "party", "ideo_fct_3", "trump_won", "existing_owner",
                                   "new_owner"),
                       label = list(
                         edu_ ~ "Education",
                         urban_type ~ "Community",
                         income_brackets_4 ~ "Income Bracket",
                         ideo_fct_3 ~ "Ideology",
                         trump_won ~ "Trump won",
                         party ~ "Party",
                         age ~ "Age",
                         race ~ "Race",
                         gender ~ "Gender",
                         existing_owner ~ "Pre-Existing Owner",
                         new_owner ~ "New Owner"
                       ),
                       missing = "no")


out <- tbl_merge(list(out1, out2, out3, out4), c("**All Respondents**", "**Pre-Existing Gun Owners**",
                                                 "**Pandemic Purchasers**", "**New Owners**")) %>%
  modify_header(
    update = list(
      label ~ "**Variable**"
    )) %>%
  modify_footnote(
    update = all_stat_cols() ~ "Reweighted"
  ) %>% 
  modify_caption("**Table A1: Sample Demographics**")

```

```{r}
gtsave(out %>% as_gt(), filename = "appendix_tables/a1.png")
```

